---
title: "Natural experiment Iraq"
author: "Dai Yamao"
date: "`r format(Sys.time(), '%Y-%m-%d')`"
output: html_document
---

```{r setup, include=FALSE}
rm(list=ls())
knitr::opts_chunk$set(echo = TRUE, fig.width = 10, fig.height = 5)
require(memisc)
require(stringi)
library(stargazer)
library(coefplot)
library(interplot)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(stargazer)
library(stringr)
library(tidyverse)
library(ggplot2)
library(psych)
#library(margins)
library(dplyr)
library(modelsummary)
library(makedummies)
library(ggeffects)
library(estimatr)
library(ggpubr)
library(BalanceR)
library(scales)

theme_set(theme_sjplot())

# Iraq
dat_iraq <- read_csv("data/Iraq.csv")
dat_iraq <- dat_iraq[-1,]

dat_iraq$RecordedDate <- as.Date(stri_match_first_regex(dat_iraq$RecordedDate, "\\d{4}-\\d{2}-\\d{2}"))

dat_iraq$ceasefire <- NA
dat_iraq$ceasefire[dat_iraq$RecordedDate <= as.Date("2025-01-16")] <- "Pre"
dat_iraq$ceasefire[as.Date("2025-01-17") <= dat_iraq$RecordedDate] <- "Post"
dat_iraq$ceasefire <- factor(dat_iraq$ceasefire)

dat_iraq <- dat_iraq |>
  mutate(Gender = D1,
         Age = D2,
         Education = D3,
         Income = D5,
         Fatah_support = Q4_3)

dat_iraq <- dat_iraq |> 
  mutate(Shia_party_supporter = (Q4_1 + Q4_2 + Q4_3)/3) |>
  mutate(Shia_party_supporter = ifelse(Shia_party_supporter >= 7, "Shia  Party Supporter", "Others"))
dat_iraq$Shia_party_supporter <- as.numeric(dat_iraq$Shia_party_supporter)

# Lebanon
dat_lebanon <- read_csv("data/Lebanon.csv")
dat_lebanon <- dat_lebanon[-1,]

dat_lebanon$RecordedDate <- as.Date(stri_match_first_regex(dat_lebanon$RecordedDate, "\\d{4}-\\d{2}-\\d{2}"))

dat_lebanon$ceasefire <- NA
dat_lebanon$ceasefire[dat_lebanon$RecordedDate <= as.Date("2025-01-16")] <- "Pre"
dat_lebanon$ceasefire[as.Date("2025-01-17") <= dat_lebanon$RecordedDate] <- "Post"
dat_lebanon$ceasefire <- factor(dat_lebanon$ceasefire)

dat_lebanon <- dat_lebanon |>
  mutate(Gender = D1,
         Age = D2,
         Education = D3,
         Income = D5,
         Hizbualla_support = Q4_3)

# Palestine
dat_palestine <- read_csv("data/Palestine.csv")
dat_palestine <- dat_palestine[-1,]

dat_palestine$RecordedDate <- as.Date(stri_match_first_regex(dat_palestine$RecordedDate, "\\d{4}-\\d{2}-\\d{2}"))

dat_palestine$ceasefire <- NA
dat_palestine$ceasefire[dat_palestine$RecordedDate <= as.Date("2025-01-16")] <- "Pre"
dat_palestine$ceasefire[as.Date("2025-01-17") <= dat_palestine$RecordedDate] <- "Post"
dat_palestine$ceasefire <- factor(dat_palestine$ceasefire)

dat_palestine <- dat_palestine |>
  mutate(Gender = D1,
         Age = D2,
         Education = D3,
         Income = D5,
         Hamas_support = Q4_2)

# Yemen
dat_yemen <- read_csv("data/Yemen.csv")
dat_yemen <- dat_yemen[-1,]

dat_yemen$RecordedDate <- as.Date(stri_match_first_regex(dat_yemen$RecordedDate, "\\d{4}-\\d{2}-\\d{2}"))

dat_yemen$ceasefire <- NA
dat_yemen$ceasefire[dat_yemen$RecordedDate <= as.Date("2025-01-16")] <- "Pre"
dat_yemen$ceasefire[as.Date("2025-01-17") <= dat_yemen$RecordedDate] <- "Post"
dat_yemen$ceasefire <- factor(dat_yemen$ceasefire)

dat_yemen <- dat_yemen |>
  mutate(Gender = D1,
         Age = D2,
         Education = D3,
         Income = D5,
         Huthi_support = Q4_1)

```


# timeseries
```{r}
Sys.setlocale("LC_TIME", "en_US.UTF-8")

# IQ
dailycount_iq <- dat_iraq |>
  group_by(RecordedDate) |>
  summarise(Count = n()) |>
  arrange(RecordedDate)
dailycount_iq$RecordedDate <- as.Date(dailycount_iq$RecordedDate)

line_date <- as.Date("2025-01-16")

p1 <- ggplot(data = dailycount_iq, aes(x = RecordedDate, y = Count)) +
  geom_bar(stat = "identity", fill = "darkgrey") +
  geom_vline(xintercept = as.numeric(line_date), color = "blue", linetype = "dashed", size = 1) +
  annotate("text", x = line_date, y = max(dailycount_iq$Count) * 1.05, 
           label = "Ceasefire", color = "red", angle = 0, vjust = 0) +
  scale_x_date(date_breaks = "5 days", date_labels = "%b %d") +
  labs(title = "Sample per day in Iraq",
       x = "Date",
       y = "Number of samples") 
p1

# LB
dailycount_lb <- dat_lebanon |>
  group_by(RecordedDate) |>
  summarise(Count = n()) |>
  arrange(RecordedDate)
dailycount_lb$RecordedDate <- as.Date(dailycount_lb$RecordedDate)

p2 <- ggplot(data = dailycount_lb, aes(x = RecordedDate, y = Count)) +
  geom_bar(stat = "identity", fill = "darkgrey") +
  geom_vline(xintercept = as.numeric(line_date), color = "blue", linetype = "dashed", size = 1) +
  annotate("text", x = line_date, y = max(dailycount_lb$Count) * 1.05, 
           label = "Ceasefire", color = "red", angle = 0, vjust = 0) +
  scale_x_date(date_breaks = "5 days", date_labels = "%b %d") +
  labs(title = "Sample per day in Lebanon",
       x = "Date",
       y = "Number of samples") 
p2

# PL
dailycount_pl <- dat_palestine |>
  group_by(RecordedDate) |>
  summarise(Count = n()) |>
  arrange(RecordedDate)
dailycount_pl$RecordedDate <- as.Date(dailycount_pl$RecordedDate)

p3 <- ggplot(data = dailycount_pl, aes(x = RecordedDate, y = Count)) +
  geom_bar(stat = "identity", fill = "darkgrey") +
  geom_vline(xintercept = as.numeric(line_date), color = "blue", linetype = "dashed", size = 1) +
  annotate("text", x = line_date, y = max(dailycount_pl$Count) * 1.05, 
           label = "Ceasefire", color = "red", angle = 0, vjust = 0) +
  scale_x_date(date_breaks = "5 days", date_labels = "%b %d") +
  labs(title = "Sample per day in Palestine",
       x = "Date",
       y = "Number of samples") 
p3

# Ym
dailycount_ym <- dat_yemen |>
  group_by(RecordedDate) |>
  summarise(Count = n()) |>
  arrange(RecordedDate)
dailycount_ym$RecordedDate <- as.Date(dailycount_ym$RecordedDate)

p4 <- ggplot(data = dailycount_ym, aes(x = RecordedDate, y = Count)) +
  geom_bar(stat = "identity", fill = "darkgrey") +
  geom_vline(xintercept = as.numeric(line_date), color = "blue", linetype = "dashed", size = 1) +
  annotate("text", x = line_date, y = max(dailycount_ym$Count) * 1.05, 
           label = "Ceasefire", color = "red", angle = 0, vjust = 0) +
  scale_x_date(date_breaks = "5 days", date_labels = "%b %d") +
  labs(title = "Sample per day in Yemen",
       x = "Date",
       y = "Number of samples") 
p4

# plot
p5 <- ggarrange(p3, 
                p2 + theme(axis.title.y = element_blank()),
                p1, 
                p4 + theme(axis.title.y = element_blank()),
                nrow = 2, ncol = 2, align = "hv")
p5

ggsave(filename = "result/sample_per_day.png",
       plot = p5,
       width = 15,
       height = 7,
       units = "in", 
       dpi = 600,
       bg = "white")
```

# balance check: Iraq
```{r}
balance_iq <- BalanceR(data = dat_iraq, 
                    group = ceasefire,
                    cov = c(Gender, Age, Education, Income))
print(balance_iq, digits = 3)
summary(balance_iq, digits = 5)
plot(balance_iq, point.size = 5, text.size = 18)

balance_IQ <- plot(balance_iq, point.size = 5, text.size = 18) +
  labs(title = "Iraq")
ggsave(filename = "result/balance_IQ.png",
       plot = balance_IQ,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)
```

# balance check: Lebanon
```{r}
balance_lb <- BalanceR(data = dat_lebanon, 
                    group = ceasefire,
                    cov = c(Gender, Age, Education, Income))
print(balance_lb, digits = 3)
summary(balance_lb, digits = 5)
plot(balance_lb, point.size = 5, text.size = 18)

balance_LB <- plot(balance_lb, point.size = 5, text.size = 18) +
  labs(title = "Lebanon")
ggsave(filename = "result/balance_LB.png",
       plot = balance_LB,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)
```

# balance check: Palestine
```{r}
balance_pl <- BalanceR(data = dat_palestine, 
                    group = ceasefire,
                    cov = c(Gender, Age, Education, Income))
print(balance_pl, digits = 3)
summary(balance_pl, digits = 5)
plot(balance_pl, point.size = 5, text.size = 18)

balance_PL <- plot(balance_pl, point.size = 5, text.size = 18) +
  labs(title = "Palestine")
ggsave(filename = "result/balance_PL.png",
       plot = balance_PL,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)
```

# balance check: Yemen
```{r}
balance_ym <- BalanceR(data = dat_yemen, 
                    group = ceasefire,
                    cov = c(Gender, Age, Education, Income))
print(balance_ym, digits = 3)
summary(balance_ym, digits = 5)
plot(balance_ym, point.size = 5, text.size = 18)

balance_YM <- plot(balance_ym, point.size = 5, text.size = 18) +
  labs(title = "Yemen")
ggsave(filename = "result/balance_YM.png",
       plot = balance_YM,
       width = 7,
       height = 4,
       units = "in", 
       dpi = 600)
```

# plot
```{r}
p10 <- ggarrange(balance_PL,
                 balance_LB + theme(axis.title.y = element_blank()),
                 balance_IQ,
                 balance_YM + theme(axis.title.y = element_blank()),
                 nrow = 2, ncol = 2, align = "hv")
p10
ggsave(filename = "result/balance_natural_experiment.png",
       plot = p10,
       width = 15,
       height = 8,
       units = "in", 
       dpi = 600)
```

# regression: Iraq
```{r}
reg01 <- lm(formula = Fatah_support ~ ceasefire, 
            data = dat_iraq)
summary(reg01)

reg02 <- lm(formula = Fatah_support ~ ceasefire +
              Gender + Age + Education + Income, 
            data = dat_iraq)
summary(reg02)

p_iraq <- plot_model(reg01, type = "pred", terms = c("ceasefire"),
           axis.labels = "ceasefire", title = "Iraq")
```

# regression: Lebanon
```{r}
reg10 <- lm(formula = Hizbualla_support ~ ceasefire, 
            data = dat_lebanon)
summary(reg10)

reg11 <- lm(formula = Hizbualla_support ~ ceasefire +
              Gender + Age + Education + Income, 
            data = dat_lebanon)
summary(reg11)

p_lebanon <- plot_model(reg10, type = "pred", terms = c("ceasefire"),
           axis.labels = "ceasefire", title = "Lebanon")
```

# regression: Palestine
```{r}
reg20 <- lm(formula = Hamas_support ~ ceasefire, 
            data = dat_palestine)
summary(reg20)

reg21 <- lm(formula = Hamas_support ~ ceasefire +
              Gender + Age + Education + Income, 
            data = dat_palestine)
summary(reg21)

p_palestine <- plot_model(reg20, type = "pred", terms = c("ceasefire"),
           axis.labels = "ceasefire", title = "Palestine")
```

# regression: Yemen
```{r}
reg30 <- lm(formula = Huthi_support ~ ceasefire, 
            data = dat_yemen)
summary(reg30)

reg31 <- lm(formula = Huthi_support ~ ceasefire +
              Gender + Age + Education + Income, 
            data = dat_yemen)
summary(reg31)

p_yemen <- plot_model(reg30, type = "pred", terms = c("ceasefire"),
           axis.labels = "ceasefire", title = "Yemen")
```

# plot
```{r}
dq_result <- ggarrange(p_lebanon, p_palestine,
                       p_iraq, p_yemen,
                       nrow = 2, ncol = 2, align = "hv")
dq_result

ggsave(filename = "result/direct_question_result.pdf",
       plot = dq_result,
       width = 8,
       height = 5,
       units = "in", 
       dpi = 600)
```